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In the present communication we consider the one-dimensional (ID) isotopically disordered lattice 
with the harmonic potential. Our analytical method is adequate for any ID lattice where potential 
energy can be presented as the quadratic form U = \ JV . q(i) Uij q(j), where q(i) - coordinate or 
velocity of i-th particle. There are derived the closed system of equations for the temporal behavior 
of the correlation functions. The final expressions allow to calculate the kinetics and dynamics 
of the system - energy, temperature profile, thermal conduction and others. There is developed 
the method for the calculation of the evolution of the eigenvalues (frequencies) and eigenvectors 
(relaxation times) to their stationary values. Exact results are obtained for times ~ 10 14 . The 
methods are suggested allowing to extend the range of the relaxation times upto ~ 10 28 . The 
spectrum of relaxation times reaches it constant value starting from the number of particles N in 
the lattice N > 300. Thermal conductance n has the non-monotonic character: for the number 
of particles N < 300 k increases as k ~ 2.4 lg N, reaches the maximum value equal ~ 4.0 at 
N ~ 300 and then slowly decreases upto V = 700. The stationary state is unique and satisfies the 
Gibbs distribution. An excellent agreement between numerical simulations and analytical results is 
obtained where possible. 

PACS numbers: 05.45.Pq, 05.50+q, 05.60.Cd 



At present, more attention is focused on statistical and transport properties of low-dimensional 
systems^. Besides obvious achievements there exists few unresolved fundamental problems, e.g. if 
there exists the unique stationary state? If this state exists, then what is the time of relaxation 
to stationary state ? Are some values self-averaged, i.e. if the thermodynamical limit really exists 
for these values? What are the sufficient conditions for the observation of finite value of thermal 
conductance in low-dimensional systems? 

The diverging value of thermal conductance k oc N a with 0.17 < a < 0.5 was obtained for 
certain model systems^^^l. Theoretical estimations for a give values 2/5&& or 1/3^ depending 
on the chosen model. 

There was obtained the final value of thermal conductivit}i£ii2*ii*i£, and, moreover, numerical 
simulations predict even the "phase transition" from normal to diverging thermal conductivity 
when temperature and/or parameters vary. However some of these results were criticized^. 

Usually it is supposed that non-integrability is the necessary but not enough condition for the 
final value of thermal conduction in ID systems^. There was demonstrated^ that any ID system 
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with the acoustical branch of excitation should have an infinite thermal conductivity in the limit 
of low temperatures. This statement is supported by the theorem that the one-dimensional system 
with momentum conservation and nonzero pressure have infinite thermal conductivity^. 

Many model lattices with different potentials of interaction were investigated. Important qual- 
itative and quantitative results were obtained, though controversial results were observed for the 
same model systems. These and associated problems are thoroughly considered in the recent 
Review^. 

ID lattices are very useful prototypes for the (analytical and numerical) investigation of kinetic, 
dynamical and transport properties of more complex and practically interesting systems such as 
carbon nanotubes (see e.gi&H^^). 

The problem of the dynamics of the disordered lattices was formulated over 50 years ago by 
F.DysonSS, who gave the general formulation of the problem. 

Dyson's problem was chosen for the thorough investigation as this system meet most of the prob- 
lems more or less common in the simulation of thermal conductance in low-dimensional systems. 
Small number of parameters allows to analyze their influence on the final results, and dynamical 
and transport properties. 

We start from the quadratic hamiltonian with fixed boundary conditions (we are mainly inter- 
ested in methodology, and the role of boundary conditions and other parameters will be considered 
and presented in a separate publication.) 

N N 

H=-J2m(i)v 2 (i) --J2x{i)U ijX (j), (1) 

i=l j=l 

where in the considered case the matrix U has the form: 



U = 



(-2 1 ••• \ 

1 -2 1 ••• 

1 -2 ••• 

••• 1 -2 1 

\ ••• 1-2/ 



(2) 



Recall that our approach is valid for the arbitrary quadratic Hamiltonians of the form 

H = '£q(i)H ij q(j), (3) 

giving rise to the linear equations of motion, where q(i), q(j) - coordinates or velocities of the 
particles. 

For the modelling of the heat reservoir we use the Langevin random forces with friction, where 
in addition to dynamical equations the members [— ^(i) x(i) +£(«)] were added, and £(i) - 

the "friction" coefficient and random force, correspondingly. Lets consider first the most general 
case, when the Langevin sources acts on all particles of the lattice. Then the equation of motion 
in matrix form is 



q = 



(Aq + S), (4) 
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where time dependent vector of the state q is represented for convenience in the form 

q(t) = ( J ) = Ml), "(2), ■ • • v(N);x(l), x(2), . . . , x(N)}, (5) 

vector = {£(1), £(2), ...,£(N)} in (0} is the vector of uncorrelated random forces, where 
random forces obey the standard relations: 

h) C0"> *a)) = 2 «(ti - ta) 7(0 T(i)/m 2 (i) (6) 

Matrix j4 has the block form 

and matrix [7 namely this one that used in (JIJ and J5J. 

In equation J7J M(i,j) — Sij m(i) is the diagonal matrix of random masses of particles, and 
M _1 - inverse matrix. We have chosen well known from literature the mass distribution such 
that randomly and with the equal probability masses have values 1 or 1/2. Matrix T — 5y - 
diagonal matrix of friction coefficients in Langevin forces for the corresponding particles, matrices 
/ and 0, are, correspondingly, unit and null matrices. 

Generally speaking, matrix A falls into a category of random matrices widely used in mathe- 
matics and physics (see, e.g. recent reviewSi). 

Initial condition is chosen in the form q(t = 0) = 0, i.e. initial temperature of all particles are 
equal to zero. Then the formal solution of equation J3J has the form 

q{t)= I W(t-r)E(r)dT (8) 
Jo 

Note, that if the initial conditions are arbitrary, i.e. q(t = 0) = <fo, then there the solution of the 
inhomogeneous solution should be added to the r.h.s. of © and it has a form: 

q(t) = W{t) q + f W(t- t) S(t) dr. (9) 
Jo 

The time evolution operator W{t) satisfies the equation 

^=AW, W(t = 0)=I (10) 

Define now the 2N x 2N matrix of correlation functions (keeping in mind the definition (JSJ ) 

(q(i, t) Q(.h t)) = f f dr x dr 2 W iM (t - Tl )W 3 . n (t - r 2 ) ( n ) & (r 2 )> . (11) 
Jo Jo 

Using the expression ©, one can get: 

C ti = (q(i, t)q(j, t)) = E 27 ^) (P) Jl drW ip (r)W jp (r). (12) 
where p is the number of a particle in the lattice subjected to Langevin forces. 
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Expression for Wi p (r) describes the evolution of the hole system, when at t = the velocity of 
p-th particle is equal to unity, and displacement from the equilibrium position is zero. All other 
particles have zero displacements and velocities. 

In the general case the evolution operator Wi p (t) for the p-th particle forecasts its temporary 
behavior at the initial condition Wi P (t = 0) = 5i p : 

dWi, 



ip 



dt 

3 

Generally speaking, one can solve the system of 2 iV equations l|13fl for every of Langevin particles 
p. As it can be seen from Q12|L the contributions from Langevin sources are also independent and 
additive. It also means the linearity of results on temperatures of Langevin sources. 

The standard approach for the solving the system (|13[) consists in transformation to the problem 
for eigenvalues. Namely this approach was used in this communication. The solution will be 
searched as the series by eigenfunctions of operator A. 

Let's consider the problem for the eigenvalues of relaxing vibrational states (vibrations relax 
because of an existence of friction in Langevin forces): 

Aq k = X k q k , (14) 

where k is the number of eigen solution (vibration). The complex value X k can be represented in 
the form: 

X k = -— + iuj k , (15) 

T k 

where r k and tu k are, correspondingly, the relaxation time and eigen frequency of fc-th mode. X k 
is not purely imaginary because of the relaxation in Langevin forces, and formally because the 
matrix A is the unsymmetrical matrix and its eigenvalues in the general case are complex values. 
Let's expand Wi P in the series over eigenvalues q k : 

2N 
k=l 

Actually we have N equations for all Langevin particles with still unknown coefficients c kp , which 
can be determined from the same equation 1)16(1 at t = and the initial condition Wi P {t = 0) = 5i p : 

2N 

S v = z2 c k P q(i)k, (17) 
what in the symbolic form can be rewritten as 

B = Ax, B = S ip , x = c kpi A = q k (i). (18) 

Substituting the solution of the equation (|16fl in <|12() , one finally gets 

c » = 1^ — 2T~r~ 1^ — ^ — T\ q q fc 2 (19) 
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At t — > oo the fraction in (|19|l reduces to l/(\ kl + Afc 2 ) 

One can check, that our results ((12(1 and (|19f) satisfies the equation for the correlation function, 
which can be obtained from the Fokker-Planck equation for the distribution function P(q) with 
the initial condition P(q,t = 0) = 5(q). 

An expression 112(1 is easily generalized for the correlation functions taken at different instants 
of time: 

(q(i, t) q(j, t + A)) = ^2 2 TS r y f drW ip (r)W jp (r + A) (20) 
p=1 m \P) Jo 

The generalization of expression l|19|) for correlation functions taken at different times is also 
obvious. 

Further particular calculations will be done for the case, when the Langevin forces acts only on 
extreme 1-st and N-th particles and we'll present some results for this case. We also consider for 
the definiteness that 7(1) = -f(N) — 1. The solution of this problem we find as the solution of the 
problem for eigenvalues (|14J) and the system of linear equations (11711 . 

Note, that the standard matrix calculus with double precision (at the diagonalization of matrix 
Q) give the relaxation times only of the order r ~ 10 14 . Really, as the result we get the complex 
value, and if the imaginary part (frequencies) ~ 1, then it is impossible to "catch" the real part 
1/t < 10~ 14 . This fact becomes obvious if one returns back to the expression (15} ■ There stands 
a value of the type e luJt in the nominator, and if co is defined with the precision ~ 10~ 14 , then it 
is not possible to use times greater then t > 10 14 . Nevertheless we can use t = 00, because of the 
damping the phase terms becomes to be zero. 

Namely these long-living states correspond to highly localized states. It follows from the well 
known fact, that there are localized states with eigen frequencies ~ 1 and very small damping in 
long disordered chains. Analogous localized states were discovered by P. AndersonSS in the diffusion 
problem (see also review 2 ^ ). 

Actually, in double summation 1)19(1 the most dangerous is the case, when in denominator A^ = 
A£ 2 , what is the very small number for the localized vibrational modes ~ —2 t p . 

But it turns out that the order of computations can be doubled. Namely, from the equation 
Aq = Xq, rewritten in components, it follows 

m(p) x k (p) A 2 . = x k (p- 1) - 2x k (p) + x k (p + 1) - 7(p) X k x k (p). (21) 

From the equation 1(211) one can get (multiplying on complex conjugated, summing and dividing 
real and imaginary parts) expression (valid only at 7(1) = 'y(N) = 1): 

1_ = M1)| 2 + MAQ| 2 

Tk 2 YstLl m i \ X i\ 2 

Thus there stands a value of the order of unity in denominator. If eigenvector is known with the 
precision of the order of 10~ 14 , then the relaxation times can be calculated with the precision 
~ 10~ 28 , i.e. t ~ 10 28 . It means that using the standard PC one can investigate the kinetics of 
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FIG. 1: The dependence of normalized distribution function for relaxation times for different number of 
particles in the lattice. Solid curve - N = 100, dashed line - N = 200. Points of different forms correspond 
to N from N = 320 to N = 700. 



isotopically disordered chains for the times much greater then it is available at any direct numerical 
simulation. 

Thus, the limiting (by precision) stages of computations are the matrix operations Q and l|18|) . 
Below we shall demonstrate some of the most essential results. 

The spectrum of relaxation times depending on the chain length N at N > 100 is shown in 
Fig. 1. As we mainly interested in the long-time kinetics, then the calculation of the relaxation 
times we've performed starting from lgr > 7.5 (up to this value the spectrum is irregular and it 
is approximately the end point of relaxation times for the regular harmonic chain). One can see, 
that the spectrum depends on N for N < 300. Starting from N = 300 the spectrum takes the 
constant value, as in the thermodynamical limit this spectrum should be the constant value. 

The "tails" of relaxation times were approximated by several functional dependencies /(lg(r)), 
and the mean of logarithm of relaxation time in all cases was calculated according to 

L°t lgr f(lgr) dlgr , v 

It turned out that at this approach and the different choices of approximating functions (lg r) 17. 
Naturally, we supposed that the relaxation law is valid for larger times, therefor the upper limit in 
(|23[1 is chosen to be the infinity. 

For N > 300 the probability that lgr > 7, 5 is large then 1/3. i.e. area of the tail is large then 
1/3 from the total area under the curve of the distribution function. 

The largest time have vibrational modes, localized in the center of the chain. It is connected with 
the fact that the corresponding vibrational wave functions have negligibly small amplitudes at the 
Langevin sources (ends of the chain), and these sources too long "drive" these modes. Naturally, 
that the maximum relaxation times T max should grow with the growth of N. But the law of growth 
of T m ax with N is still unknown. 
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The temperature profiles along the chain for differen boundary temperatures are shown in Fig. 2 
for N = 70, where our results are exact, for arbitrary times. The profiles are shown for r = 
4, 6, 10, 14, oo time units. 

As was shown above, because of extremely long life time of vibrational modes, localized in 
the center of the lattice and having the relaxation times greater than > 10 28 for N > 100, the 
temperature in the center of the lattice is calculated not very exactly. For the lattice of any size 
we can guarantee only the results for ~ 30 — 40 particles from every end. It also means that 
the numerical simulations (e.g., temperature distribution along the chain) for large lattices are 
incorrect as the states localized in the center of the lattice have no enough time to relax to the 
stationary state during the period of time of the numerical simulation. 

Because of the thermal flow in the stationary state is a constant value along the length of all 
the lattice, then it can be calculated in the vicinity of ends, where, as was stated above, all our 
calculations are exact. It is obvious that the system may not reach its stationary state. The 
specific coefficient of thermal conductance k as dependent on the chain length is shown in Fig. 3. k 
depends only on the temperature difference on the chain ends and doesn't depend on the absolute 
temperature ((T(l) + T(N))/2), what is obvious from (|12fl . Really, if the Langevin sources acts at 
extreme atoms (1-st and iV-th), then the expression (|12fl can be rewritten 

2T(1) f l 2T(N) f* 

Ct ^m^T) I W ^ W ^ dT + ^W) I W lN (r)W ]N (r)dr. (24) 

Let's change the temperature at the ends in the same manner: T(l) — > T(l) + T, T(N) — > 
T(N) + T. Then the change in correlation functions l|24[) has the form of the correlator for equal 
temperatures at the ends of the lattice: 

Oij = ^ I W ll {T)W jl { T )dT + ^ T f W tN (r)W jN (T)dT, (25) 

m l JO m N JO 
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FIG. 3: The dependence of the coefficient of thermal conductance on the chain length. The linear fitting 
of the initial part of the curve corresponds to k ~ 2.4 lg N. 

and it is obvious, that any changes in the thermal flow and other additives, depending on the 
temperature difference become zero. 

As it is seen from Fig. 3, for not too long lattices, k grows as k ~ 2,4 lg TV upto N ~ 300, and 
then slowly decreases. (Calculations for N > 1000 we couldn't perform because of computational 
limitations). 

Note, that all phonon states are localized for the ID latticaSi^. It was proven for the electro- 
conductance problem 24 that the electroconductivity tends to zero as the absolute temperature also 
tends to zero. As to our knowledge, there is no such a statement for the thermal conductance of 
disordered ID systems. But in the weak coupling approximation, i.e. at 7 <C 1 it was supposed 
that for the rigid boundary conditions n ~ 1/%/iV, i.e. tends to zerrA (Recall the strange result 
that at free boundary conditions thermal conduction diverges as n ~ yN~) . 

It is probable that the slow decreasing of k for N > 300 obsereved in Fig. 3 is the manifestation 
of the tendency of k to zero when N — ► 00. We couldn't check this fact because of computational 
limitations. 

The thermal conduction was calculated according to k = — J AT, where AT was defined as the 
temperature difference on the lattice ends (it is impossible to determine the temperature in the 
vicinity of a lattice center, because of an existence of localized modes with extremely long time 
of life) . Therefor we restricted ourself by such a simple expression which is valid by the order of 
magnitude. 

In Conclusion we say that the solved problem partially continues the classical works of Libowtz 
with co-workers2&^i on the dynamics of ID disordered systems. We investigated the kinetic, dy- 
namical and transport properties of ID isotopically disordered harmonic lattice. The accuracy of 
our computations for all values, except the relaxation times, is limited by values of the order of 
10 -14 . We've succeeded to calculate the relaxation times with the precision 10 -28 . Our precision 
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of solution is determined by standard matrix operations performed on PC with double precision. 
Relaxation times calculated by us many orderers of magnitude exceed any value obtained in di- 
rect numerical simulations. We've got non-monotonic dependence of thermal conduction which 
is dependent on N with maximum at N ~ 300. The stationary state is unique but times of its 
achievement grow with the increase in chain length. Because of additivity of all results on tem- 
perature (jl9|) . the detailed calculations can be done at single chosen temperature at the chain 
ends. 

This work was supported by the RFBR (02-02-16205, 04-03-32119 , 04-02-17306) and the Pro- 
gram of Presidium of Russian Academy of Sciences " Fundamental problems of physics and chem- 
istry of nanosized systems and nanomaterials" . 
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